# Create "FloodInsTakeup" variable

# Setup -----
library(tidyverse) 
library(tidycensus)



########################################################################################
# Download & load NFIP flood insurance policies data 
# SOURCE: https://www.fema.gov/about/openfema/data-sets, "OpenFEMA Dataset: FIMA NFIP Redacted Policies"
# Download and load data as "NFIP_Redacted_Policies"


NFIP_Redacted_Policies$policyEffectiveDate   <- as.Date(NFIP_Redacted_Policies$policyEffectiveDate)   # convert to date format
NFIP_Redacted_Policies$policyTerminationDate <- as.Date(NFIP_Redacted_Policies$policyTerminationDate) # convert to date format


# Extract number of NFIP policy-in-force as of 2017-12-31 by census tract -----------
NFIP_Redacted_Policies_FL = NFIP_Redacted_Policies %>% 
  filter(propertyState=="FL") %>%                                # only keep property in FL 
  filter(policyEffectiveDate<=as.Date("2017-12-31") &            # Effective date <= 2017-12-31
         policyTerminationDate>as.Date("2017-12-31")) %>%        # Termination date > 2017-12-31
  group_by(censusTract) %>%
  summarise(NFIP_unit = sum(policyCount))                           # Number of policy-in-force as of 2017-12-31 by census tract 



# Extract the 2017 estimate of the number of housing units by census tract from ACS 5-year dataset
tract2017 <- get_acs(geography = "tract", 
                  variables = c(HousingUnit = "B25001_001"), 
                  state = "FL", 
                  year=2017,
                  survey = "acs5")

tract2017 = tract2017 %>% select(GEOID, HousingUnit2017 = estimate)



# Combine and calculate a tract's NFIP flood insurance takeup rate as of 2017-12-31

NFIPTakeup = tract2017 %>%
  left_join(NFIP_Redacted_Policies_FL, by=c("GEOID"="censusTract")) %>%
  filter(HousingUnit2017>0) %>%
  mutate(
    FloodInsTakeup = ifelse(is.na(NFIP_unit), 0, NFIP_unit/HousingUnit2017)
    )



save(NFIPTakeup, file="/02_Data/NFIPTakeup.Rdata")


